Bouw LHM topologie
Inlezen LKM25 links¶
Lees in kolom nodefrom en nodeto: lkm25\Schematisatie\KRWVerkenner\shapes\LKM25_Links.shp
- vinden LSW eindpunten
- Per eindpunt LKM-links naar lopen tot node_to start met LSM
- Zoek SOBEK
- Wanneer LSW ID gelijk is aan nodefrom: volg het netwerk totdat nodeto start met LSM.
- Check nu of er een Sobek lateral op minder dan 250m. ligt. Zoja verbindt deze LSM met deze SOBEK lateral.
Zoniet loop het netwerk verder af en check nogmaals of het volgende eindpunt wel op 250m. van een lateral ligt. Let op: Je mag alleen Sobek laterals gebruiken die in hetzelfde district liggen als de LSW.
from config import DATA_DIR, LKM25_DIR, MOZART_DIR, LHM_DIR, LSW_DIR, LSM_KOPPELING_DIR, load_src
import geopandas as gpd
load_src()
from lhm.read import read_lsm_lhm, read_lsw_routing, read_dw_keys
from lhm.lsm import snap_to_waterbodies
from lhm.lsw import lsw_end_nodes, lsw_network
Inlezen LSWs¶
Vanuit de geleverde lsws.shp en lswrouting.dik maken we een netwerk van LSWs.
lsw_routing_dik = MOZART_DIR / r"mozartin/lswrouting.dik"
lsw_routing_df = read_lsw_routing(lsw_routing_dik)
lsw_gdf = gpd.read_file(LSW_DIR / "lsws.shp")
lsw_gdf = lsw_gdf.dissolve(by="LSWFINAL").reset_index()
lsw_links_gdf, lsw_nodes_gdf = lsw_network(lsw_gdf, lsw_routing_df)
Vinden niet-gekoppelde LSWs¶
We gaan ervan uit dat LSWs die afwateren op andere LSWs níet tevens afwateren op het DM. We gaan vanaf dit punt LSWs die niet gekoppeld zijn aan andere LSWs koppelen aan het DM. Hieronder tonen we deze LSWs op de kaart
lsw_end_nodes_gdf = lsw_end_nodes(lsw_links_gdf, lsw_gdf)
print(f"we moeten {len(lsw_end_nodes_gdf)} koppelen aan DM-knopen")
lsw_end_nodes_gdf.explore()
we moeten 2309 koppelen aan DM-knopen
Koppelen via DW-keys¶
Vanuit dwkeys.txtkunnen we per district kijken aan welke DM knoop wordt gekoppeld. Wanneer een district slechts koppelt aan 1 DM-knoop, dan betekent dit dat alle LSWs in dit district gekoppeld moeten worden aan deze knoop.
We doen dit en tonen de overgebleven nog niet gekoppelde LSWs.
dw_key_file = LHM_DIR / r"dm/txtfiles_git/dwkeys.txt"
dw_keys_df = read_dw_keys(dw_key_file)
lsw_dm_links = []
for dw, df in dw_keys_df.groupby(by=["oid"]):
dm_nodes = df.loc[df.kty == "d"].nid.unique()
if len(dm_nodes) == 1:
df = lsw_end_nodes_gdf[lsw_end_nodes_gdf["DWRN"] == dw[0]]
for row in df.itertuples():
lsw_dm_links += [(row.LSWFINAL, dm_nodes[0])]
print(f"{len(lsw_dm_links)} LSWs gekoppeld aan DM-knopen")
lsw_end_nodes_gdf = lsw_end_nodes_gdf.loc[~lsw_end_nodes_gdf.LSWFINAL.isin([i[0] for i in lsw_dm_links])]
print(f"We moeten nog {len(lsw_end_nodes_gdf)} koppelen aan DM-knopen")
lsw_end_nodes_gdf.explore()
1784 LSWs gekoppeld aan DM-knopen We moeten nog 525 koppelen aan DM-knopen
lsm3_locations_csv = LSM_KOPPELING_DIR / "LSM3_locations.csv"
knoop_district_csv = LSM_KOPPELING_DIR / "LSM3_DMKnoopDistrict_childs.csv"
lsm_lhm_gdf = read_lsm_lhm(lsm3_locations_csv, knoop_district_csv)
lsm_lhm_gdf.explore()
Koppelen LSW via LSM-lateralen aan DM¶
Voor de nog overgebleven LSWs zoeken we de DM-knoop via LSM-lateralen. We nemen deze stappen:
- We kijken of er LSM-lateralen van het juist district liggen in het gebied van de LSW
- We kijken aan welke unieke DM-knopen deze lateralen zijn verbonden. We verbinden de LSW met al deze districten
lkm_waterlichamen_shp = LKM25_DIR / "KRW-waterlichamen_SGBP3.shp"
lkm_waterlichamen_gdf = gpd.read_file(lkm_waterlichamen_shp)
lsm_lhm_snapped_gdf = snap_to_waterbodies(lsm_lhm_gdf, lkm_waterlichamen_gdf, offset=250)
lkm_links_shp = LKM25_DIR / "LKM25_Links.shp"
lkm_links_gdf = gpd.read_file(lkm_links_shp)
ToDo: Deze code netjes verwerken in een functie
from shapely.geometry import Point
from lhm.dm import get_dm_nodes
if lsw_gdf.index.name !="LSWFINAL":
lsw_gdf.set_index("LSWFINAL", inplace=True)
init = len(lsw_dm_links)
def report_progress(iteration, total, interval=10):
percent = (iteration / total) * 100
progress = int(percent / interval)
line = '[' + '#' * progress + ' ' * ((100 / interval) - progress) + ']'
print(f'{line} {percent:.1f}% completed', end='\r')
def find_links(row, lsw_id): # algemene functie om links te vinden mbv een LKM link
lsw_dm_links = []
polygon = lsw_gdf.at[lsw_id, "geometry"]
df = lsm_lhm_gdf[lsm_lhm_gdf.within(polygon)]# we kijken of er lateralen in het gebied van de LSW liggen
if not df.empty:
lsw_dm_links +=[(lsw_id, i) for i in df.DM.unique()]
if row.nodeto.startswith("LSM"): # de link moet starten met LSM
point = Point(row.geometry.bounds[2:]) # de Point van de nodeto
waterlichaam_gdf = lkm_waterlichamen_gdf[lkm_waterlichamen_gdf.intersects(point)] # het waterlichaam waar punt op snapt
if not waterlichaam_gdf.empty: # er moet wél een waterlichaam gevonden worden
waterlichaam = waterlichaam_gdf.iloc[0] # dan pakken we het eerste waterlichaam (als het goed is is de lengte altijd 1)
df = lsm_lhm_snapped_gdf[lsm_lhm_snapped_gdf.within(waterlichaam.geometry)] # we zoeken naar lsm_lhm laterale knopen binnen het waterlichaam
df = lsm_lhm_snapped_gdf[lsm_lhm_snapped_gdf.within(waterlichaam.geometry)]
if not df.empty: # we hebben gevonden!
lsw_dm_links = [(lsw_id, i) for i in df.DM.unique()] # hier maken we links van de lsw_id naar de unieke DM-knopen
lsw_dm_links = [i for i in lsw_dm_links if i[1] in dm_nodes_stringified]
return lsw_dm_links
all_links = []
for idx, row in enumerate(lsw_end_nodes_gdf.itertuples()):
report_progress(idx, len(lsw_end_nodes_gdf))
lsw_id = row.LSWFINAL
district = row.DWRN
dm_nodes = get_dm_nodes(dw_keys_df,district)
dm_nodes_stringified = [str(i) for i in dm_nodes]
lkm_links_iter = lkm_links_gdf[lkm_links_gdf.nodefrom == str(lsw_id)].itertuples()
for row in lkm_links_iter:
links = []
links = find_links(row, lsw_id)
if (not links):
counter = 0
while (not links) and (counter < 10) :
df = lkm_links_gdf[lkm_links_gdf.nodefrom == row.nodeto]
if not df.empty:
row = df.iloc[0]
links = find_links(row, lsw_id)
counter += 1
all_links += [i for i in links if i not in all_links]
all_links
lsw_dm_links += all_links
print(f"{len(all_links)} LSWs koppelingen gemaakt aan DM-knopen")
lsw_end_nodes_gdf = lsw_end_nodes_gdf.loc[~lsw_end_nodes_gdf.LSWFINAL.isin([i[0] for i in all_links])]
print(f"We moeten {len(lsw_end_nodes_gdf)} koppelen aan DM-knopen")
lsw_end_nodes_gdf.explore()
457 LSWs koppelingen gemaakt aan DM-knopen We moeten 193 koppelen aan DM-knopen
Koppelen LSW aan dichtsbijzijnde DM-knoop¶
Als uiterste mogelijkheid koppelen we LSW-knopen aan de meest dichtbijzijnde DM-knoop die voorkomt in de DM-keys
Locatie DM knopen inlezen
# Define the shapefile paths
DM_nodes_shp = LHM_DIR / r"dm/data_dvc/DM_nodes.shp"
DM_links_shp = LHM_DIR / r"dm/data_dvc/DM_links.shp"
# Read each shapefile separately
DM_nodes_gdf = gpd.read_file(DM_nodes_shp)
DM_links_gdf = gpd.read_file(DM_links_shp)
DM_nodes_gdf.columns = [i.strip() for i in DM_nodes_gdf.columns]
DM_links_gdf.columns = [i.strip() for i in DM_links_gdf.columns]
# Calculate the closest distance LSW endpoints and DM
all_links = []
for row in lsw_end_nodes_gdf.itertuples():
lsw_id = row.LSWFINAL
district = row.DWRN
dm_nodes = get_dm_nodes(dw_keys_df,district)
dm_nodes = [str(i) for i in dm_nodes]
dm_node = DM_nodes_gdf.at[DM_nodes_gdf.loc[DM_nodes_gdf.ID.isin(dm_nodes)].distance(row.geometry).sort_values().index[0], "ID"]
all_links += [(lsw_id, int(dm_node))]
lsw_dm_links += all_links
print(f"{len(all_links)} LSWs koppelingen gemaakt aan DM-knopen")
lsw_end_nodes_gdf = lsw_end_nodes_gdf.loc[~lsw_end_nodes_gdf.LSWFINAL.isin([i[0] for i in all_links])]
print(f"We moeten {len(lsw_end_nodes_gdf)} koppelen aan DM-knopen")
193 LSWs koppelingen gemaakt aan DM-knopen We moeten 0 koppelen aan DM-knopen
from shapely.geometry import LineString
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
lsw_dm_links_gdf = gpd.GeoDataFrame(lsw_dm_links, columns = ["node_from", "node_to"], geometry = gpd.GeoSeries(), crs=28992)
def make_line_string(row):
report_progress(row._name, len(lsw_dm_links_gdf))
point_from = lsw_nodes_gdf.loc[row.node_from]
point_to = DM_nodes_gdf.loc[DM_nodes_gdf.ID == str(row.node_to)].geometry
return LineString([[point_from.x, point_from.y], [point_to.x, point_to.y]])
lsw_dm_links_gdf.loc[:, "geometry"] = lsw_dm_links_gdf.apply((lambda x: make_line_string(x)), axis=1)
lsw_dm_links_gdf["node_from"] = lsw_dm_links_gdf["node_from"].astype(str)
lsw_dm_links_gdf["node_to"] = lsw_dm_links_gdf["node_to"].astype(str)
[######### ] 100.0% completed
Resultaat¶
In de kaart hieronder zie je het resultaat
import folium
#m = lsw_gdf.explore(name="lsws", color="lightgrey")
m = lsw_links_gdf.explore(name="lsw links", color="blue")
m = lsw_dm_links_gdf.explore(m=m, name="lsw-dm links", color="green")
m = DM_links_gdf.explore(m=m, name="dm links", color="red")
folium.LayerControl().add_to(m)
m